Computing an unbiased inverse

In a previous entry, Reweighting is biased, I had discussed the problem of a finite-sample size bias in reweighting. In this entry I will investigate a possible solution.

But first, I had asked if this was known in the literature. A comment by Lee Warren pointed to this article, Population size bias in descendant-weighted diffusion quantum Monte Carlo simulations, which discusses this very issue (in a slightly different context, but the root issue is the same). In that context the problem is even more severe since the straightforward solution (increase N) carries a heavier computational cost.

One of the crucial issues is finding a good estimator for the inverse (1/D). The naive estimator (simply do the division) is biased. So the question is: are there any other ways of performing division that could be used? Yes, in digital arithmetic there are a couple of algorithms for division using only multiplication and addition - Newton-Raphson and Goldschmidt (see the Wikipedia entry on digital division)

In general, to have an unbiased estimator, each sample of the random variable must enter the estimator linearly (alternately, we can say each sample should only be used once and then discarded)

The Newton-Raphson iteration for the inverse (1/D) is x_i+1 = x_i(2-D x_i). To make this unbiased, in each iteration the value of D needs to be an independent sample, and the values of x_i from previous iterations also need to be computed from independent samples. One can think of the algorithm forming a tree, with each iteration branching to two children (each child is one of the x_i computed by the previous iteration), until reaching the leaves representing the initial iteration (x₀ = T₁ + T₂ D)

I tested this algorithm, and it appears to give unbiased results. There is a catch, though. The result has a large variance. Using the average of all the samples with the naive estimator yields a bias that is much smaller than the variance of the unbiased estimator (hence the mean squared error is a more appropriate metric for the combination of variance and bias)

Open questions/issues:

• Is there some modification of the algorithm that will produce a smaller MSE than the naive estimator? Or even an MSE that isn't too much larger?


• How many iterations are necessary - how to analyze the convergence of this stochastic Newton-Raphson iteration?


• Try out the Goldschmidt algorithm.